Surface-hopping dynamics and decoherence with quantum equilibrium structure 
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In open quantum systems decoherence occurs through interaction of a quantum subsystem with its 
environment. The computation of expectation values requires a knowledge of the quantum dynamics 
of operators and sampling from initial states of the density matrix describing the subsystem and bath. 
We consider situations where the quantum evolution can be approximated by quantum-classical 
Liouville dynamics and examine the circumstances under which the evolution can be reduced to 
surface-hopping dynamics, where the evolution consists of trajectory segments evolving exclusively 
on single adiabatic surfaces, with probabilistic hops between these surfaces. The justification for 
the reduction depends on the validity of a Markovian approximation on a bath averaged memory 
kernel that accounts for quantum coherence in the system. We show that such a reduction is often 
possible when initial sampling is from either the quantum or classical bath initial distributions. If 
the average is taken only over the quantum dispersion that broadens the classical distribution, then 
such a reduction is not always possible. 



I. INTRODUCTION 

The exact quantum dynamics of processes such as pro- 
ton transfer in chemical or biological systems vibra- 
tional energy relaxation in condensed phases,^ or ultra- 
fast dynamics of photoexcited molecular systems^ can- 
not be computed directly, due to the exponential scal- 
ing of computational time with the number of degrees 
of freedom of the system. Consequently, a number of 
different approximate quantum dynamical schemes have 
been used to circumvent this difficulty. Surface-hopping 
schemes j2i2jI£ where classical degrees of freedom evolve 
on single adiabatic surfaces and make probabilistic hops 
from one surface to another, are the most commonly used 
such methods. Generally, one would expect the system 
to evolve into a coherent state as a result of the coupling 
between adiabatic states. The justification for the restric- 
tion of the evolution to single adiabatic surfaces is based 
on the fact that interaction with the environment rapidly 
destroys coherence resulting in the collapse of the system 
onto a single adiabatic surface. Several prescriptions have 
been put forth to capture this physicsi 9 ' 10 ! 11 ] 12 ! 13 ! 14 ? 15 
Mixed quantum-classical Liouville dynam- 
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can be used to model quantum 
dynamics coupled to a classically evolving environment, 
and naturally accounts for nonadiabatic transitions 
and quantum coherence. Recently we showed how this 
theory could be recast into the form of a master equa- 
tion, so that the dynamics involves classical trajectory 
segments on single adiabatic surfaces, interrupted by 
hops between these surfaces - dynamics akin to that 
in surface-hopping schemes* 2 ^ In this formulation, the 
transition probabilities are determined by coherent evo- 
lution trajectory segments where two adiabatic surfaces 
are coupled. Decoherence, which is again attributed to 
interaction with the environment, is accounted for by 
averaging the coherent evolution segments over an initial 
distribution of the bath degrees of freedom. 



This analysis was applied to the computation of the 
quantum reactive flux correlation function and rate con- 
stant. The reaction rate coefficient was computed in an 
approximation where the reaction coordinate and bath 
equilibrium distributions were taken to be classical. Here 
we study the effects that result from sampling from quan- 
tum equilibrium distributions. This allows us to study 
a number of issues related to quantum initial states, 
decoherence and the validity of simple surface-hopping 
schemes. 

The outline of this article is as follows. In the next sec- 
tion we present a brief overview of the reduction of the 
quantum-classical Liouville equation to a master equa- 
tion description, along with a discussion of the computa- 
tion of off-diagonal matrix elements of a bath averaged 
operator in terms of the evolution of the diagonal ele- 
ments. In Sec. Mil we consider a model for a reactive bar- 
rier crossing process in order to illustrate the effects of 
decoherence on the reaction rate using classical and quan- 
tum initial sampling. We compare the results of mas- 
ter equation and quantum-classical Liouville dynamics. 
The quantum initial distribution can be written in terms 
of a classical distribution broadened by quantum disper- 
sion. In Sec. IIVI we investigate whether averages over 
the quantum dispersion are sufficient to justify a Marko- 
vian description of the dynamics. Finally, we present our 
conclusions. 



II. THEORY 

The average value of a quantum mechanical operator 
is given by 



A{t) = TvAp{t) = TrA(t)p(0), 



(1) 



where A(t) is some observable whose time evolution is 
given by the Heisenberg equation of motion and p(0) is 
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the initial quantum density matrix. Depending on the 
problem of interest, it may be convenient to represent 
a subset of the degrees of freedom in some basis. In 
this study we find it convenient to use an adiabatic ba- 
sis but any other basis may be chosen. The remaining 
degrees of freedom are Wigner transformed,— providing 
a phase-space like representation of this part of the sys- 
tem. The choice of which degrees of freedom are Wigner 
transformed is based on a consideration of the problem of 
interest. Thus, an equivalent exact quantum expression 
for the average value in Eq. ([T|) is 



A{t) = Tr' / dXA w (X,t)p w (X), 



dXA^ (X,t) Pw a (X), 



(2) 



where the adiabatic states and energies are given by 
the solutions to the eigenvalue problem, hw{R)\ot; R) = 
E a (R)\a;R), with H w = hw{R) + 537- In these ex- 
pressions, X denotes phase space variables (R,P), Tr' 
denotes the partial trace over the remaining quantum de- 
grees of freedom and the subscript W denotes the partial 
Wigner transform over the selected degrees of freedom. 
The evaluation of this expression involves sampling from 
an initial quantum distribution followed by the evolution 
of the dynamical variable. 

Here, as in previous studie o 26 i 27 we assume that the full 
quantum evolution can be replaced by mixed quantum- 
classical Liouville evolution given b y 20 ' 28 



— A aa ' 



Mr 



{X, t) = Y] iC aa -,w A w (*> *)• ( 3 ) 



The quantum-classical Liouville super-operator, C 



aa',0/3' 



(4) 



The two operators in the first term on the right hand 
side are the frequency corresponding to the energy gap, 
u aa > = (E a — E a ')/h, and classical Liouville operator 

iL aa > = Ttf ■ -§5 + \i F w + F w) ' -§pi respectively. The 
Hellmann-Feynman force for state a is F w . The last 
term, 3 a a' fip' is the operator responsible for nonadia- 
batic transitions and corresponding momentum adjust- 
ments which ensure energy conservation. The complete 
definitions of these operators can be found in fiefs. [2(| 
and (H- 

We define the subsystem as those degrees of freedom 
that are of primary interest, while the bath comprises the 
remainder of the system. For example, in the study of 
proton transfer—, the subsystem was chosen to be the 
protonic degree of freedom, q, plus the solvent polar- 
ization, while the bath consisted of the remaining sol- 
vent and molecular complex degrees of freedom. Let- 
ting X = (Rq,P ) and Xb = (Rb,Pb) be the vari- 
ables in the subsystem and bath, respectively, that are 



Wigner transformed, the Hamiltonian can be written as, 
Hw — Ho + -fffe(o) > where the subsystem and bath-plus- 
interaction Hamiltonians are given by, 



H = 



H 



6(0) 



P 

2m 

EL 

2M 
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2M 



v s (q,Ro), 



Vbo(Rb] Rq)- 



(5) 
(6) 



For this Hamiltonian, the adiabatic states depend upon 
the subsystem coordinates, Rq, while the adiabatic ener- 
gies, E a (R), depend upon the full set of coordinates. 

The subsystem density matrix, p s (Xo), and the bath 
density, p c (Xb] Xq), conditional upon the subsystem con- 
figuration, i?o, are defined as 



p s [X ) = j dX b p w (X ,X b ), (7) 
p c (X b ;X ) ee Pw (X ,Xb)p: 1 (X ). (8) 

We can express the initial quantum density in the adia- 
batic basis as, 

P&(x ,x b ) = Y / pf(x b ;X )pP/(x ). (9) 



Using this definition, we may rewrite Eq. @ as, 

A$j = E/ dX (A w (X,t))r'p^' a (X ). (10) 

aa' 

This expression is now composed of an average of the dy- 
namical variable over the conditional bath distribution, 

(A w (X,t))r' =E / dX b A^(X,t)p^'(X b ;Xo), 

followed by an average over the subsystem density. The 
computation of this quantity thus consists of sampling 
from the initial subsystem quantum distribution, and 
evolution of the bath averaged dynamical variable. In 
the Schrodinger picture, one would expect that averaging 
the density matrix over the bath distribution would cause 
the coherences, the off-diagonal elements of the subsys- 
tem density matrix, to decay rapidly. This idea served as 
motivation to derive a master equation, which provided 
a simple surface-hopping description of the dynamics. 



A. Quantum-classical master equation 

The quantum-classical master equation was derived by 
writing Eq. ([3]) as a coupled set of evolution equations for 
the diagonal and off-diagonal elements. By formally solv- 
ing this set of coupled differential equations, one obtains 
a generalized master equation for the evolution of the 
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diagonal elements of the operator 



24 



-A«(X 1 t) = iL a A a d {X,t) 



dt' 



^M^(X,t')A^X^ t „X b!t/ ,t- t') 




where the memory function, M^a(X,t') is defined as, 



M${X,t) = 2Rc [W a p(t',O)}D al3 {X o )D a0 (X Oa p tt ,). 



(13) 

Here D a /3(Xo) — (Pq/Mq) ■ d a p(Ro), the nonadiabatic 
coupling matrix is d a /3(R ) = (a; R \V Rq \/3; R ), W a (3 is 
a phase factor, and the subscripts and superscripts on 
the memory function label the indices on the first and 
second D functions respectively. For the remainder of 
this article, superscripts and subscripts that appear on 
variables denote evolution along a particular adiabatic 
surface, (aa'), and the appearance of both a subscript 
and superscript on a variable indicates the action of two 
consecutive nonadiabatic transitions. A bar appearing 
over Xq, denotes a shift in the momentum coming from 
the J operator. The shift is applied to the momentum 
coordinate, Pq, such that^ 



and ja^p is a momentum shift operator^ Since Eqs. (JSj) 
and (|10[) are equivalent, computation of the average value 
(12) using either Eq. {2} and the master equation after the 
lift to the full phase space or Eq. (|10[) with the non- 
Mar kovian subsystem equation will yield the same result. 
Note that the subscripts on the second memory term in 
Eq. (115j) are the same. This term arises from the memory 
function corresponding to (M£™(X, t))b- Trajectories ac- 
counted for by this term jump to the mean surface and 
then return to their original surface. Thus, the net effect 
is no jump, but a phase factor is introduced. 

The master equation (fT5|) provides a trajectory de- 
scription that prescribes evolution on single adiabatic 
surfaces interspersed with quantum transitions between 
them. The probability of a transition is obtained entirely 
from the coherent evolution. Furthermore, decoherence 
is accounted for by the decay of the off-diagonal evolution 
segments resulting from averaging over the bath distri- 
bution. 



Pq = Po + d a p(sgn(P ' d a p) 



x V (P • d al3 ) 2 + AE a pM - (P • d af3 ) , 



B. Off-diagonal elements 

The expression for the average value in Eq. (|10[) in- 
cludes a sum over all of the matrix elements of the observ- 
able. Here we show that the evolution of the off-diagonal 
matrix elements can be expressed in terms the diagonal 
matrix elements so that the master equation can be used 
(14) to approximately compute these contributions. 

The off-diagonal part of the quantum-classical Liou- 
villc evolution equation ([3]) is given by 



and Xq = (Po,Po). Given that the operator is initially 
diagonal, the evolution equation (| 1 2[) is formally equiva- 
lent to the mixed quantum-classical Liouvillc equation, 20 
only now it is organized in two parts: evolution on single 
adiabatic surfaces and the memory kernel containing all 
of the off-diagonal evolution. 

A generalized master equation was derived by project- 
ing Eq. (|12p onto the subsystem space and making a 
Markovian approximation on the bath averaged mem- 
ory kernel that accounts for the coherence in the system. 
The bath average provides a mechanism for decoherence 
and leads to decay of the memory kernel. The resulting 
equation is still non-Markovian in character as a result 
of the projection of the adiabatic evolution onto the sub- 
system. The final master equation is obtained by lifting 
the equation back into the full phase space to recover a 
fully Markovian master equation description: 



dt 



A a d (X,t)=iL a A a d {X,t) 



(15) 



- J2 m a p(X o )j a ^ A d (X, t) - m aa (X )A a d (X, t) 





where 



m a 0(X o ) 



dt'(M^(X,t')) b , (16) 



-A (X,t) = i£°A (X,t)+i£°' d A d (X,t), (17) 
at 

where the subscripts d and o indicate diagonal and off- 
diagonal matrix elements respectively. Formally solving 
this equation gives 

A {X,t) = e lCOt A o (X,0) 

+ [ dt l e iC ° { - t 'kC°' d A d (X,t-t l ). (18) 



If we assume that the observable is initially diagonal, 
then the first term vanishes. Applying the definition of 
C given in Eq. Q, the second Liouville operator in the 
above expression reduces to i£ o d = —J o d . Substituting 
into Eq. (fill) we nnd 



A aa '(X,t) = [ dt' ]T U o aa ,^,(X,t')J^, A^X,t-t% 
Jo w',0 



_ l e iC°(X)t^ 



(19) 



. aa, is the off- 



where U^ a , )/30l (X,t) 

diagonal propagator responsible for evolution in off- 
diagonal space. For a two level quantum system, this 
propagator is given exactly by 

K a >,00>{t) - W aa >{t,Qy L ™ l{ - x)t 5 a p5 a ,p, , (20) 
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for a, a', (3, (3' = 1,2 and a ^ a',f3 ^ This can be 
verified by using the Dyson formula and noting that, for a 
two level system, the only off-diagonal propagator matrix 
elements that contribute to the dynamics are U° 2 12 — 
U$* 2i ■ This definition is generally applicable for weak 
nonadiabatic coupling. Acting on the observable with 
this propagator and using the definition 20 - 24 of J the 
evolution equation becomes 

A aa '{X 1 t)= [ dt'MZ a '(X ,t')A aa '(X,t,t'), (21) 
Jo 

where we have introduced the off-diagonal memory func- 
tion, 

M™'{X ,t) = W aa ,(t)D aa ,(X$), (22) 



and 



A aa '(X,t,t') = (A2'(X$,X b , t ,,t-t') (23) 

- A d(^0,t' > X b ,f,t ~ t') 



The expression given in Eq. (|2ip prescribes evolution 
along the mean surface for a time t' followed by a jump to 
a diagonal surface. Thus, once again, this memory func- 
tion is made up entirely of coherent evolution segments 
which carry phase factors. 

The off-diagonal memory function differs from that in 
Eq. (I13|) by a single momentum-jump operator. In the 
diagonal form, evolution starts on an adiabatic surface, 
jumps to the mean surface and then jumps to another 
adiabatic surface; thus, two jump operators are involved. 
In contrast, off-diagonal evolution starts on a mean sur- 
face and, therefore, requires only one jump to get to an 
adiabatic surface. 

By taking the bath average of Eq. (|2"Tj) and factoring 
the average on the right hand side, we obtain an evolution 
equation for the bath averaged observable in the reduced 
subsystem phase space: 



(A aa (X,t)) b 



dt'(Mr'(X ,t')) b (A aa ') b (t,t'). 



(24) 



The arguments which are used to justify this approxima- 
tion have been given in the Appendix of Ref . [24] . The 
bath average of the oscillating phase factor in the mem- 
ory function once again leads to decay. If this decay is 
rapid, it provides justification for making a Markovian 
approximation, 



(Mr' (x ,t')) b 



dt'(M« a '(X ,t')) b )5(t') 



2m* a '(Xo)5(t'), 



(25) 



so that the expression for the evolution of the off-diagonal 
matrix elements of an observable is given by 



Lifting the equation back to the full phase space we have 
A aa \x,t)= (27) 
mr'(Xo) (j aa ,(X )A%'(X,t)-j a , a A%(X,t)) . 

The dynamics of the diagonal terms is given by Eq. (|15p , 
and the momentum shift operators, j Qa < , defined in 
Ref. [24J, have been pulled out of the observable. We 
have obtained an expression for the evolution of the off- 
diagonal matrix elements of an observable entirely in 
terms of diagonal evolution. This expression prescribes 
a transition from the mean to each adiabatic surface at 
time t = with probability m"" (Ao), followed by diag- 
onal evolution. The resulting trajectories are the same 
as those given by the diagonal master equation only the 
momentum is shifted at time t = 0. 

Since upward jumps from the mean surface require en- 
ergy input to the subsystem, one must check that there 
is sufficient energy in the bath to provide the subsys- 
tem with the required energy. Transitions that satisfy 
AE a pMo/ (Pq ■ dap) 2 > 1 are not allowed. In the case 
where this requirement is satisfied, the off-diagonal evo- 
lution is only a function of the evolution along the lower 
adiabatic surface. 



III. REACTION RATE AND QUANTUM 
SAMPLING 

The master equation theory was used to compute the 
chemical reaction rate for a process A ^± B assuming that 
the bath and subsystem reaction coordinate were treated 
classically^ Here, we consider the quantum mechanical 
expression for the rate coefficient that involves averag- 
ing over the quantum equilibrium structure. This allows 
us to investigate the implications of quantum equilibrium 
sampling on decoherence and the quantum-classical mas- 
ter equation approximation to the reaction rate. 

The forward rate constant is given by2£ 



kAB (t) 



1 



(2 - S a ' a ) 



(28) 



dARc 



N% a (X t t)W% 



where W% a (A, ^ ) is the spectral density function that 
depends on the quantum equilibrium structure, and 
Ng a (A, t) is the time evolved matrix element of the 
number operator for the product state B. To calculate 
the rate, one samples initial configurations from quan- 
tum equilibrium distributions, and then computes the 
evolution of the number operator for product state B. 

The reaction model we consider here is the same as 
used in our previous study^ The subsystem consists of 
a two-level quantum system bilinearly coupled to a quar- 
tic oscillator and the bath consists of v — 1 = 300 har- 
monic oscillators bilinearly coupled to the non-linear os- 
cillator but not directly to the two-level quantum system. 



(A aa (X,t)) b = m™ (X )(A aa (X,t)) b . 



(26) 
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The Hamiltonian and system parameters are the same as 
those used in prior studies ; 24 ! 31 except the potential pa- 
rameters are a = 0.25, and b = 1.0, (3 — 1/ksT is taken 
to be j3 = 1.0 and (3 = 0.5 corresponding to low and high 
temperatures, respectively, and the nonadiabatic cou- 
pling strength is 70 = 1.25. The reaction coordinate is Rq 
and the product species operator N^ a (Ro) = 9(Ro)S aa ' 
is initially diagonal in the adiabatic basis. Here 9(Rq) is 
the Heaviside function. 

For this reaction model, the barrier region of the po- 
tential is locally harmonic along the reaction coordinate 
Rq. Consequently, one can construct an approximate ex- 
pression for the spectral density function that has the fac- 
torized form, WA{X,ihj3/2) ss pa(X )pI(Xi,; Ro), where 
Pa(Xq) is the subsystem spectral density function and 
p%{Xb] Rq) is the Wigner transform of the canonical equi- 
librium density matrix for the bath in the field of the Rq 
coordinate. 30 Inserting this form of the spectral density 
function into the rate coefficient expression, the second 
line of Eq. (f2"5|) may be replaced by 



dXRe 



N r'(x,t)w2' a (x, 1 -^) 



dX Re 



{ Nr') b (Xo,t)pi a (Xo, l -f-) 



(29) 



where the angle brackets indicate an average over the 
conditional equilibrium distribution, p^(Xi,',Ro). The 
computation of the rate constant now involves the bath 
averaged observable discussed in Sec. HH allowing us to 
make use of the formalism discussed above. 



A. Simulation Results 

The evolution of the observable A^" (X, t) was com- 
puted using both quantum-classical Liouville dynam- 
ics and master equation dynamics. The simulation 
of quantum-classical Liouville dynamics was carried 
out using the sequential short-time propagation algo- 
rithm 3 ^ in conjunction with the momentum-jump ap- 
proximation^ 8 .^ and a bound on the observable . 29 ! 34 
The initial positions and momenta of the quartic os- 
cillator were Monte Carlo sampled from the harmonic 
part of the quantum equilibrium density distribution 
/c/4 "(Xq, ih(3/2)21 and the bath coordinates were sampled 
from the quantum equilibrium density for harmonic os- 
cillators with appropriate frequencies. Simulation of the 
master equation consisted of two stages: calculation of 
the memory kernel, followed by the sequential short-time 
propagation algorithm restricted to single adiabatic sur- 
faces^ Quantum transitions were Monte Carlo sampled 
using the appropriate memory function for the transition 
probabilities. 

Computation of the off-diagonal memory function in- 
volves a procedure similar to that of the diagonal mem- 
ory function^ A bath configuration is sampled from the 
quantum equilibrium density for a given subsystem co- 



ordinate, A"o, and then a trajectory is evolved adiabati- 
cally along the mean surface for a prescribed time t. The 
memory function is calculated by numerically integrating 
Eq. (f2"2"|) along the trajectory and then averaging over an 
ensemble of initial bath configurations. This was done 
for a range of Xq values resulting in a memory function 
surface. The surface generated by integrating this quan- 
tity is plotted as a function of Xq in Fig. Q] for the low 
temperature case. The surface for the high temperature 
case is very similar in structure; however, there are minor 
numerical differences that one would expect to see due 
to weaker quantum effects at higher temperatures. 




FIG. 1: Plot of m (Ro, P ) versus Ro and P for (3 = 1.0. 

Once the memory functions are known, we may cal- 
culate the evolution of the observable as discussed else- 
where.— The off-diagonal contribution may be calculated 
using Eq. (|27[) . The initial bath configuration and sub- 
system coordinate are sampled from the off-diagonal el- 
ements of the full quantum equilibrium structure as be- 
fore; however, now a nonadiabatic transition from the 
mean surface to each adiabatic surface is made immedi- 
ately. If there is insufficient energy for an upward transi- 
tion, then there is no population transfer to the excited 
state for that realization. Otherwise, the trajectory con- 
tributes to the rate, weighted by m 12 (Xo). The evolution 
is then computed using the sequential short time propa- 
gation algorithm to simulate Eq. (fT5")) for each diagonal 
matrix element. 

A comparison of the reactive flux correlation function 
computed using quantum-classical Liouville dynamics, 
master equation dynamics, and adiabatic dynamics can 
be seen in Fig. [2] The lowering of the rate due to nonadi- 
abatic effects is notable in both cases. The master equa- 
tion reproduces the short time build-up of the quantum 
correlation function in both cases in spite of the fact that 
the short time build-up occurs on the same time scale as 
the decay of the bath averaged memory kernel. The zero 
initial value of the correlation function is a consequence 
of the quantum mechanical treatment of the reaction co- 
ordinate. As expected, the error bars associated with 
the master equation calculation are smaller than those 
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FIG. 2: Forward rate constant kAB(t) as a function of time. In 
both figures, the upper curve is the adiabatic rate, the middle 
curve is the quantum master equation result and the lowest 
curve is the mixed quantum-classical Liouville result. In both 
panels the microscopic relaxation time t m i C is approximately 
3.4 units. Ra)1ff = 1.0. RTjJI/3 = 0.5. 



from the quantum classical Liouville simulation since the 
master equation calculation does not involve oscillating 
phase factors. 

The reaction rate constant is given by the plateau val- 
ues of the correlation function. The full rate computed 
with the master equation is close to but larger than the 
quantum-classical Liouville calculation. From Fig. [3] one 
sees that the discrepancy is due to the off-diagonal and 
excited state contributions to the rate. 

The quantum-classical Liouville and master equation 
ground state diagonal contributions to the rate agree very 
well. There are discrepancies between the two methods 
in the excited state and mean surface contributions. For 
the excited state contribution, the discrepancy can be 
understood in terms of recrossing of the barrier region. 
In the master equation description, trajectories starting 
in the excited state, jump down to the ground state, ac- 
companied by an increase of momentum, and rapidly sta- 
bilize. In contrast, for quantum-classical Liouville evolu- 
tion, trajectories starting in the excited state first make 
a transition to the mean surface; subsequent evolution 
results in recrossings on this surface, followed by a tran- 



FIG. 3: Ground state, excited state and off-diagonal contri- 
butions to the forward rate constant fcAs(i) as a function of 
time. The top pair of curves is the ground state contribu- 
tion, the middle pair is the excited state contribution and the 
bottom set is the off-diagonal contribution. In each pair, the 
smoother curve is the result from the master equation sim- 
ulation while the noisier curve is from the quant urn- classical 
Liouville calculation. In both cases it is clear that the dis- 
crepancy arises from the disagreement in the off-diagonal and 
excited state contributions. |(a)|/3 = 1.0. |(b)|/? = 0.5. 



sition to the ground state where it stabilizes. This results 
in increased recrossing of the barrier region, lowering the 
plateau value of the correlation function. Although the 
ground state contribution is subject to the same quali- 
tative dynamics, upward transitions are associated with 
a reduction in the momentum and thus stay confined to 
the barrier region in both descriptions. As a result, there 
is a smaller difference between the recrossing corrections 
for the ground state contribution in the two descriptions. 

The deviations between the off-diagonal terms are 
largely a consequence of the approximations made in 
Sec. Ill Bl In particular, the Markov approximation gives 
rise to a non-zero initial value, which ultimately results 
in a higher plateau value. Similar effects have been seen 
in the context of quantum master equations i 35 i 36 Short 
time transients typically relax before the slow subsys- 
tem processes take place, and the neglect of this behav- 
ior leads to a shift in the initial conditions which affects 
the longer time value of computed quantities. The rele- 
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vant time scale to consider in relation to this effect is the 
decay rate of the memory function or the decoherence 
time, Tdecoh ■ The probability distribution of decoherence 
times was calculated from the first zero crossing of the 
bath average of the memory function for an ensemble of 
Xq values. The results of this calculation can be seen 
in Fig. [H The mean decoherence time is slightly longer 
for the low temperature case with Tdecoh ~ 0.992, while 
the high temperature value is Tdecoh ~ 0.738. The dif- 
ference can be attributed to the form of the distribution 
of decoherence times, which exhibits a long tail in the 
lower temperature case. The time scale of the decay to 



time that characterizes the relaxation of the correlation 
function to the plateau value, ranges between 10 and 4 for 
our simulation results. In order to provide some insight 
into these values, we remark that our reaction model pro- 
vides a simple representation of proton transfer in solu- 
tion. 51 In a more realistic model of proton transfer, t lnic 
was found to be approximately 300 fs. 29 Using the scaling 
reported here, our decoherence time would correspond to 
tens of femtoseconds, a value that is reasonable for con- 
densed phase systems. A precise assignment of real time 
units depends on how the model parameters are mapped 
into those of physical systems. 




FIG. 4: Distribution of the decoherence times for the bath 
averaged memory function at (a) high and (b) low tempera- 
tures. 

the plateau value in the off-diagonal contribution to the 
correlation function is of the same order of magnitude as 
T decoh- Thus, one can infer that the discrepancy appear- 
ing in this contribution is due to the neglect of transient 
memory effects. 

The decoherence time as defined above characterizes 
the decay associated with the bath average of the dy- 
namics of off-diagonal density matrix elements in the adi- 
abatic basis. Recall that the memory kernel results from 
substituting the solution of the equation of motion for 
the off-diagonal density matrix elements into the evolu- 
tion equation for the diagonal elements. A number of dif- 
ferent measures of the decoherence time have been used 
in the literature. A basis-independent quantity that has 
been used to determine the decoherence time is Tr'p 2 , 
where p s is again the subsystem density matrix averaged 
over the bathi^I This quantity contains information re- 
garding population transfer in addition to the decay of 
off-diagonal density matrix elements, resulting in biexpo- 
nential decay. The time associated with the fast decay 
segment of the curve can be identified with the decoher- 
ence time^ We have computed the decoherence time for 
the high temperature parameter set from the short time 
decay of Tr'p 2 and obtained a value of Tdecoh ~ 0.80 in 
agreement with that obtained for the decay of the mem- 
ory kernel. We note here that the ratio of decoherence 
and microscopic time scales t m i C / Tdecoh, where t m i C is the 



IV. QUANTUM DISPERSION 



The above results were obtained by sampling from 
an equilibrium quantum distribution while in our ear- 
lier study of decoherence and surface-hopping dynam- 
ics 2 ^ equilibrium initial states were sampled from a clas- 
sical initial distribution. We now explore the differences 
between the two descriptions. 

The quantum mechanical treatment of the reaction co- 
ordinate has the most significant effect on the structure of 
the reactive flux correlation function* 3 ^ The correlation 
function rapidly increases from zero, (Fig. [2]), and this is 
qualitatively different from the classical case where the 
reactive flux correlation function has a non-zero initial 
value equal to the transition state theory rate. A quan- 
tum mechanical treatment of the reaction coordinate has 
a more significant effect on the reaction rate than a quan- 
tum treatment of the bath. 39 However, we are interested 
in decoherence, which arises from a bath average for fixed 
values of the reaction coordinate. Therefore, it is impor- 
tant to examine the consequences of treating the bath 
quantum mechanically or classically. 

The memory surface, which is used to compute the 
nonadiabatic transition rates, is calculated by averaging 
the memory function over an initial bath distribution. 
In Fig. [5] we compare the surfaces obtained using the 
quantum and classical bath distributions. The surface 
obtained using quantum sampling is somewhat broader 
than the classical surface. This is expected since quan- 
tum dispersion will broaden the distribution. The magni- 
tude of the variation is quite small and has a small effect 
on the transition rates. 

A different perspective on how the quantum bath 
structure affects decoherence can be obtained by noting 
that the quantum equilibrium distribution for a bath of 
harmonic oscillators can be written as a convolution of 
the classical density, pf(Xb]R ) = l/Zfe~^ Hb , with its 
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FIG. 5: Memory surface calculated using (a) the quantum 
bath equilibrium distribution, and (b) classical bath equilib- 
rium distribution. 



quantum dispersion to obtain 
d 



dt 



(A3{X',t)) g (X,t) = (iL a A<Z(X>,t)) g (X,t) (33) 



dt'{M(X', t')A d (X', t - t')) g {X, t - t'), 



where we have simplified the notation of the memory 
term analogous to that in Appendix B of Ref. [HI ■ Since 
the factored form of the bath distribution involves a con- 
volution, one cannot use the projection operator formal- 
ism that was employed in the derivation of the master 
equation. Instead, we assume that the observable varies 
slowly over the width of the quantum dispersion function, 
so that it may be taken out of the integral. Furthermore, 
if the observable is approximately constant over the inte- 
gration region then, A%(X,t) « (A°ftX',t)} g (X,t). Us- 
ing these approximations we obtain 



quantum dispersion function ; 39 ! 40 

.IXlgiXi-X^ptiX^Rol (30) 



Pb(Xb; Rq) 
g{x b ) 



— TT 



Z h -Ll 2t:(u'' - 1) 



xe 



, (31) 



where, 



ij coth Uj , and Uj 



Phuj/2, The effect of 
the quantum dispersion function, g{X' b — Xj,), a Gaus- 
sian centered at the bath phase point X b , is to broaden 
the classical density around each point in phase space. 
This function accounts for the quantum fluctuations of 
the bath. We consider whether averages over the bath 
quantum dispersion distribution of the memory kernel for 
fixed values of the reaction coordinate and bath phase 
space coordinates yield sufficiently rapid decay of the 
memory kernel to justify a Markovian approximation. 

Starting from the quantum rate expressions (|28p and 
(|29|) . we substitute the convoluted form of the bath den- 
sity to obtain, 

k AB (t) = / dX(N§ a '(X',t)) g (X,t) (32) 

xpi a (X , l -^)p c b l (X b ;R ), 

where the average value of the species variable over the 
quantum dispersion is denoted by (Ng a (X', t)) g (X, t) = 
J dX' b g(X' b - X b )N% a ' (X 1 , t). The resulting quantity de- 
pends on the subsystem degrees of freedom and bath 
phase space coordinates determined by the classical bath 
distribution. An approximate evolution equation for 
quantum dispersion average of an observable can be ob- 
tained a follows. 

Starting from the generalized master equation (|12p , we 
take the average of the entire evolution equation over the 



-(A%(X>,t)) g (X,t) 



(34) 



iL a (A5(X',t)) g (X,t)- / dt'{M(X', t')) g (X, t') 



x(A d (X',t-t')) g (X,t-t'). 

The quantum dispersion average of the memory function 
appears in the second term of this expression. In Fig. 6(a) 
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FIG. 6: Comparison between the quantum dispersion aver- 
age {MH(X',t)) g (X,t) and full bath average (M 1 1 |} b (f) of 
the memory function versus time, (a) For high (solid) and 



low (dashed) temperatures, the oscillations persist through 
the simulation time, (b) The full bath average of the diagonal 
memory function for high (solid) and low (dashed) temper- 
ature, for comparison to the quantum dispersion averaging. 
Note that the scale on the time axis is shorter and that the 
function decays quickly. 
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we show (M a P(X' , i)) g (X, t) as a function of time. Each 
curve in the plot represents a particular choice of X 
and Xb, averaged over the quantum dispersion function 
g(X' b — Xb). From the figure we see that the function does 
not decay on a short time scale and exhibits long time 
recurrences. This structure should be contrasted with 
that shown in Fig. |6(b)| where a full bath averaged was 
carried out. Comparing the classical distribution to the 
quantum dispersion function (see Fig. [7]), one sees that 
even at low temperatures the range of bath configurations 
described by g(X' b — X b ) is narrow. 



. 5 
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FIG. 7: Comparison of the classical bath density (solid line) 
plotted as a function of momentum at Ro — and the quan- 
tum dispersion for different oscillator frequencies. The widest 
quantum dispersion curve (dashes) corresponds to the cutoff 
oscillator frequency, Lu max — 3, while the narrow curve (dots) 
corresponds to an intermediate oscillator frequency, uij = 1.2. 
Low temperature (/3 = 1) results are plotted, where one would 
expect the quantum dispersion to be most significant. Al- 
though the curve corresponding to the maximum oscillator 
frequency is close to the width of the full bath density, it 
comprises a small contribution to the quantum dispersion av- 
erage over the distribution of oscillator frequencies. 

In this representation, each classical trajectory on the 
mean surface is replaced by an ensemble of trajectories 
whose initial values are determined by the quantum dis- 
persion. The simulation results for our reactive model 
indicate that each ensemble of trajectories centered at a 
given system phase space point is confined to a narrow 
tube for long times. Thus, decoherence is not observed on 
a time scale relevant to the decay of the reactive flux cor- 
relation function. However, if the average is taken over 
the full bath distribution, then decoherence is rapid, jus- 
tifying the Markovian approximation leading to master 
equation dynamics. 

There is another aspect of quantum dispersion versus 
full bath averaging that merits discussion. The bath pro- 
jection led to an equation of motion valid in the subsys- 
tem. The lift to full phase was a device to obtain a fully 
Markovian description to compute the subsystem expec- 
tation values^ In contrast, if quantum dispersion aver- 
aging could be used to justify a Markovian approximation 
on the memory kernel, one would obtain a master equa- 



FIG. 8: The memory surface for = 1.0, averaged over the 
quantum dispersion function using different bath configura- 
tions for each subsystem configuration. The surface is very 
rough indicat ing t hat the choice of initial bath strongly aff ects 
Memory surface for upward transitions, (b) 
'or downward transitions. 



the dynamics (a) 



Memory surface 



tion valid in the full phase space of the system. To study 
the nature of this dependence, we computed the quantum 
dispersion average of the memory surface as a function 
of the initial bath configuration. For each initial subsys- 
tem coordinate we sampled a bath configuration from a 
Boltzmann distribution and then performed an average 
over the quantum dispersion for that particular configu- 
ration. A new bath configuration was used for each sub- 
system coordinate. The results of this calculation, shown 
in Fig. [51 indicate that the memory surface is sensitive 
to the initial bath configuration. Consequently, it is not 
possible to reduce the dimensionality of the memory sur- 
face in this formulation and this makes the calculation of 
the transition probabilities difficult. The surface may be 
smoothed by averaging over the ensemble of bath config- 
urations, i.e., taking the full bath average. 



V. CONCLUSION 

Decoherence in open quantum systems has its origin 
in the interactions of the system with its environment.— 
This paper focused on the effects of quantum versus 
classical equilibrium sampling on decoherence and the 
computation of average values and correlation functions. 
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Both of these quantities can be written in forms that 
involve quantum initial sampling (or classical sampling 
when this approximation is appropriate) and quantum 
dynamics, which we model by quantum-classical Liou- 
ville dynamics. As in our earlier study— the decoher- 
ence problem can be cast in the form of the validity of a 
Markovian approximation on the bath averaged memory 
kernel for the evolution of the diagonal elements of the 
density or operators. 

The simulations on a simple reaction model with quan- 
tum initial sampling allowed us to assess the validity of 
simple surface-hopping dynamics and how it relates to 
decoherence. Since the reaction coordinate is treated 
quantum mechanically, the reactive flux correlation func- 
tion starts from zero and finally reaches a plateau value 
from which the rate constant can be determined. In spite 
of the fact that this build-up typically occurs on a similar 
time scale to the decoherence time, the initial structure 
of the correlation function is well described by master 
equation dynamics. 

The role of quantum dispersion on decoherence was 
also investigated. Depending on system parameters, av- 
erages over the quantum dispersion around fixed bath 
phase space coordinates may be insufficient to provide a 
justification for a Markovian approximation to the mem- 
ory kernel. Furthermore, simulations on the reaction 
model indicated that the quantum dispersion average of 
the memory kernel is a strong function of the bath co- 



ordinates and, thus, so are the transition probabilities in 
the master equation. 

An aim of this investigation was to examine the con- 
ditions under which simple surface-hopping descriptions 
of the dynamics, where the system evolves on single adi- 
abatic surfaces interrupted by quantum transitions be- 
tween such surfaces, are valid. An average of the mem- 
ory kernel over some distribution that results in decay 
on a time scale faster than the decay of the correlation 
function of interest is required in order for a Markovian 
master equation description of the dynamics to be valid. 
The decoherence of the subsystem must occur on a time 
scale that is shorter than that of slow subsystem pro- 
cesses. When these conditions are satisfied, the result- 
ing dynamics provides a useful tool to compute expecta- 
tion values since it is more stable with smaller statistical 
uncertainties. The results of this work should prove to 
be useful for understanding the dynamics of many body 
quantum systems more complex than those considered 
for our model calculations. 
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